In [1]:
from NDimInv.plot_helper import *
import dd_res as DD_RES
dd_res = DD_RES.dd_resistivity()
import scipy.stats as stats
Generate multiple spectra exhibiting various high-frequency characteristics, e.g. an unfinished slope or a finished peak.
In [86]:
def get_spectrum(peak1, peak2, strength1, strength2, omega, s, rho0):
"""
Generate a two-peak spectrum and chargeability distribution for a Debye Composition
peak1 and peak2 denote the mean of the norm-distributions for m values (in s coordinates).
"""
m = stats.norm.pdf(s, peak1, 1.5) * (strength1) + \
stats.norm.pdf(s, peak2, 1.5) * (strength2 * 2)
# we need to norm the chargeabilities because they will be scaled by rho0. Arbitrarily set the norm factor so that it is 1 for rho0 = 100
m = m / rho0
# we also want to normalize to the number of chargeabilities per frequency decade
f = 1 / (2 * np.pi * 10**s)
nr_per_decade = int(len(s) / (np.log10(np.max(f)) - np.log10(np.min(f))))
m = m / nr_per_decade
# apply conversion
pars = np.hstack((rho0, m))
pars_converted = dd_res.convert_parameters(pars)
re,mim = dd_res.forward_re_mim(omega, pars_converted, s)
return re, mim, rho0, m
# input parameters
rho0 = 100
Nd = 20
f = np.logspace(-3,4,30)
omega = 2 * np.pi * f
tau_s, s_s, f_s = dd_res.get_tau_values_for_data(f, Nd)
peaks = [
[0, -4.5, 'peak just inside'],
[0, -6, 'peak just outside'],
[0, -3, 'peak finished'],
[2, -1, 'peak really far']
]
re_list = []
mim_list = []
for nr,options in enumerate(peaks):
re,mim,rho0,m = get_spectrum(options[0], options[1], 1, 1, omega, s_s, rho0)
re_list.append(re)
mim_list.append(mim)
# save mag/pha to file
m_re = np.array(re_list)
m_mim = np.array(mim_list)
mag = np.abs(m_re - 1j * m_mim)
pha = np.arctan(-m_mim / m_re) * 1000
magpha = np.hstack((mag,pha))
np.savetxt('synthetic/data.dat', magpha)
np.savetxt('synthetic/frequencies.dat', f)
In [32]:
# load data
frequencies = np.loadtxt('frequencies.dat')
omega = 2 * np.pi * frequencies
data = np.atleast_2d(np.loadtxt('data.dat'))
mag = data[0, 0:data.shape[1]/2]
pha = data[0, data.shape[1]/2:]
cmplx = mag * np.exp(1j*pha/1000)
re = np.real(cmplx)
mim = -np.imag(cmplx)
odir = 'test_dd_orig'
# load final iteration
rho0 = np.atleast_2d(np.loadtxt(odir + '/rho0_results.dat'))[0,1]
#rho0 = 10**log10rho0
m = np.loadtxt(odir + '/m_results.dat')
mtot = np.sum(10**m)
#m = 10**m
s = np.loadtxt(odir + '/s_results.dat')
pars = np.hstack((rho0, m))
fre,fmim = dd_res.forward_re_mim(omega, pars, s)
del_mim_del_m = dd_res.del_mim_del_chargeability(omega, pars, s)
nr_m = del_mim_del_m.shape[1]
def plot_spectrum():
# plot spectrum, and results
fig = plt.figure(figsize=(15,6))
ax = fig.add_subplot(3,1,1)
ax.semilogx(frequencies, re, '.')
ax.semilogx(frequencies, fre, '-', color='r')
ax = fig.add_subplot(3,1,2)
ax.semilogx(frequencies, mim, '.')
ax.semilogx(frequencies, fmim, '-', color='r')
#plot_spectrum()
Imaginary part sensitivities for each chargeability
In [33]:
def plot_unnormalized():
fig, ax_list = plt.subplots(61, 1,figsize=(12,100))
for i in range(0,nr_m):
ax_list[i].semilogx(frequencies, del_mim_del_m[:,i])
fig.suptitle('Un-normalized')
fig.subplots_adjust(top=0.99,hspace=0.1)
#plot_unnormalized()
def plot_normalized():
max_sens = np.max(np.abs(del_mim_del_m))
fig, ax_list = plt.subplots(nr_m, 1,figsize=(12,100))
for i in range(0,nr_m):
ax_list[i].semilogx(frequencies, del_mim_del_m[:,i] / max_sens)
ax_list[i].set_ylim([0,1])
ax_list[i].set_xlabel('Frequency (Hz)')
ax_list[i].set_title('s: {0}'.format(s[i]))
fig.suptitle('Normalized')
fig.subplots_adjust(top=0.97,hspace=0.6)
#plot_normalized()
Now we try some comulative sensitivity values.
For each frequency, sum up all chargeability sensitivities:
$Covf(f) = \sum_{i=1}^{N_m} \left| \frac{\partial -Im(f)}{\partial m_i} \right|$
For each chargeability, sum up all frequency sensitivities
$Covm(m) = \sum_{i=1}^{N_f} \left| \frac{\partial -Im}{\partial m} \right|_f$
In [34]:
covf = np.abs(del_mim_del_m).sum(axis=1)
covf /= np.max(covf)
covm = np.abs(del_mim_del_m).sum(axis=0)
covm /= np.max(covm)
invcovm = np.abs(1 / del_mim_del_m).sum(axis=0)
invcovm /= np.max(invcovm)
fig, (ax1,ax2,ax3,ax4) = plt.subplots(4,1,figsize=(15,12))
ax1.semilogx(frequencies, covf)
ax1.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('Covf')
ax2.plot(s, covm)
ax2.set_xlabel(r'$s = log_{10}(\tau)$')
ax2.set_ylabel('Covm')
ax2.invert_xaxis()
ax3.plot(s, np.log10(10**m / mtot))
ax3.set_xlabel(r'$s = log_{10}(\tau)$')
ax3.set_ylabel(r'$log_{10}(m) / m_{tot}$')
ax3.invert_xaxis()
ax3.set_ylim([-2.4, -0.8])
ax4.semilogx(frequencies, mim, '.-')
ax4.set_xlabel('Frequency (Hz)')
ax4.set_ylabel(r'$-Im(\rho) (\Omega m)$')
fig.subplots_adjust(hspace=0.5)
Compute $\left[A^T A\right]^{-1}$ to get data (co-)variances
In [35]:
A = dd_res.del_mim_del_chargeability(omega, pars, s)
B = A.T.dot(A)
#B = np.matrix(B).I.dot(B.T.dot(B))
B = np.matrix(B).I
variances = diag(B)
vars_norm = variances / np.max(variances)
fig, (ax1,ax2) = plt.subplots(2,1, figsize=(15,5))
ax1.plot(s, vars_norm)
ax1.invert_xaxis()
a = ax2.hist(vars_norm, 100)
In [35]: